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Abstract 

Often,  rarefied  flows  of  interest  in  aerospace  applications  are  embedded  within  largely  continuum  flow  fields. 
Neither  continuum  nor  kinetic  methods  provide  both  physically  accurate  and  numerically  efficient  techniques 
to  simulate  the  entire  flow  field.  Instead,  multi-scale  methods  can  provide  the  capability  of  achieving  the 
efficiency  of  continuum  methods  in  regions  where  the  degree  of  collisional  nonequilibrium  is  small,  while 
maintaining  the  physical  accuracy  of  kinetic  methods  in  rarefied  portions  of  the  flow.  This  work  begins  with 
an  outline  of  typical  aerospace  flows  that  may  require  a  multi-scale  analysis  in  order  to  make  accurate  and 
efficient  predictions.  It  then  provides  an  overview  of  the  research  performed  in  developing  hybrid  methods 
for  partially  rarefied  gas  flows.  Finally,  some  results  derived  from  current  state  of  the  art  hybrid  codes  are 
presented  and  emerging  developments  are  described. 


1  Introduction 


Typically,  the  Knudsen  number,  shown  in  Eq.  1,  which  is  the  ratio  of  the  mean  free  path,  A,  to  a  characteristic 
length  scale,  L,  is  used  to  determine  the  degree  of  collisional  nonequilibrium  in  a  gas  flow. 

Kn=^  (1) 

When  Kn  is  large,  very  few  molecular  collisions  occur,  and  the  entire  flow  is  considered  to  be  in  collisional 
nonequilibrium,  or  rarefied.  For  rarefied  flows,  the  direct  simulation  Monte  Carlo  (DSMC)  method  [1]  can 
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provide  an  efficient  and  physically  accurate  prediction,  however  the  numerical  expense  of  DSMC  increases  as  the 
Knudsen  number  decreases.  If  Kn  <C  1  ,  the  flow  can  be  considered  approximately  in  collisional  equilibrium,  or 
continuum,  and  a  numerical  solution  of  the  Navier-Stokes  equations  with  modern  computational  fluid  dynamics 
(CFD)  can  be  found  efficiently  with  very  little  loss  in  physical  accuracy.  However,  the  physical  accuracy  of  the 
continuum,  Navier-Stokes  equations  degrades  as  the  Knudsen  number  increases. 

For  many  cases,  the  characteristic  length  scale  and/or  mean  free  path  may  vary  significantly  throughout  the 
flow  field,  especially  when  the  flow  is  in  the  transitional  regime  between  fully  rarefied  or  continuum.  In  these 
multi-scale  flows,  some  regions  may  have  large  characteristic  length  scales  compared  to  the  local  mean  free  path 
such  that  they  are  entirely  within  the  continuum  domain  and  efficient  prediction  with  kinetic  methods,  such  as 
DSMC,  is  not  possible.  These  same  flows  may  contain  regions  where  the  characteristic  length  scales  are  very 
small  compared  to  the  local  mean  free  path,  such  that  the  physical  accuracy  of  continuum  methods  breaks  down. 
Although  there  are  techniques  that  extend  the  continuum  formulation  by  use  of  slip  boundary  conditions,  they 
often  do  not  fill  the  entire  gap  of  applicability  between  continuum  and  rarefied  simulation  techniques.  Therefore, 
application  of  either  continuum  or  rarefied  simulation  methods  to  the  entire  domain  of  flow  fields  in  the  transition 
regime  can  not  maintain  both  physical  accuracy  and  numerical  efficiency. 

Alternatively,  a  multi-scale,  hybrid  method  can  be  applied  to  these  flows.  Such  a  method  uses  a  continuum 
solver  in  regions  that  are  very  near  collisional  equilibrium,  while  using  rarefied  gas  simulation  methods  in  regions 
that  exhibit  continuum  breakdown.  In  this  work,  we  use  the  definition  of  a  hybrid  particle-continuum  method 
proposed  by  Wilmoth  et  al.  [2]  as  an  algorithm  that  performs  two-way  information  exchange  between  particle 
and  continuum  simulation  methods. 

The  proceeding  section  will  outline  examples  of  multi-scale  flows  to  which  hybrid  particle-continuum  methods 
can  be  applied.  Next,  Sec.  3  outlines  the  research  that  has  been  performed  specific  to  developing  hybrid  particle- 
continuum  methods.  Section  4  gives  an  overview  of  the  published  hybrid  particle-continuum  algorithms  and 
a  selection  of  results  demonstrated  from  each  proposed  method.  Finally,  Sec.  5  provides  a  summary  of  the 
paper  and  an  outline  of  some  future  challenges  that  still  exist  to  extend  the  current  state  of  the  art  in  hybrid 
particle-continuum  simulation  techniques. 


2  Examples  of  Multi-Scale  Flows  in  Aerospace  Applications 


Flows  in  the  transitional  regime  have  small  length  scales  due  to  the  geometry  or  flow  structure.  For  example, 
the  geometric  size  of  MEMS  devices  is  sufficiently  small  that  the  molecular  nature  of  the  gas  is  important  and 
a  kinetic  treatment  may  be  required  in  some  regions  of  the  flow  to  maintain  physical  accuracy.  Conversely, 
hypersonic  vehicles  flying  in  the  upper  atmosphere  experience  flow  with  a  large  mean  free  path  due  to  the 
relatively  low  atmospheric  density.  In  addition,  these  hypersonic  vehicles  create  high  gradient  flow  features,  such 
as  strong  shocks,  strong  expansions,  and  thin  boundary  layers  that  exhibit  very  small  flow  length  scales.  A  hybrid 
particle-continuum  method  is  well  suited  to  simulate  internal  or  external  flows  that  exhibit  a  large  variation  of 
characteristic  flow  and/or  collision  length  scales.  This  subsection  will  outline  some  typical  flows  that  can  exhibit 
these  flow  features  and  may  be  suited  for  application  of  a  hybrid  particle-continuum  method. 
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2.1  Micro-Scale  Flows 

Recently,  an  increase  in  the  development  of  micro-  and  nano-electromechanical  systems  (MEMS/NEMS)  has 
occurred  due  to  the  low  costs  associated  with  batch-manufacturing,  low  energy  requirements,  and  small  size. 
Many  applications  of  MEMS/NEMS  devices,  which  include  micro-turbines  [3,  4],  micro-sensors  for  chemical  con¬ 
centrations  or  gas  flow  properties  [5,  6,  7,  8],  and  spacecraft  micro-propulsion  [9],  involve  gas  flows.  However,  the 
general  understanding  of  micro-scale  gas  flows  lags  behind  the  fabrication  and  application  of  these  MEMS /NEMS 
devices.  Experiments  have  shown  that  the  fluid  mechanics  of  micro-scale  gas  flows  are  not  the  same  as  those 
experienced  in  similar,  larger-scaled  devices.  For  example,  the  pressure  distribution  in  a  long  micro-channel  has 
been  observed  to  be  nonlinear  [10],  and  the  measured  flow  rate  in  a  micro-channel  was  higher  than  predicted  using 
the  Navier-Stokes  equations  [11].  Although  the  use  of  slip  boundary  conditions  within  a  Navier-Stokes  solver  has 
provided  some  improvement  over  application  of  the  no-slip  boundary  condition,  this  can  only  provide  a  small 
extension  of  application  of  the  Navier-Stokes  equations  in  the  near  continuum  regime.  Instead,  application  of  a 
hybrid  method  that  uses  a  physically  accurate  particle  method  in  all  regions  where  the  Navier-Stokes  equations 
are  invalid  (rather  than  just  at  the  surface)  can  provide  reliable  predictions  throughout  the  rarefied-continuum 
flow  regime. 


2.2  Hypersonic,  Blunt-Body  Flows 

At  moderately  high  altitudes,  some  regions  of  the  flow  about  blunted,  hypersonic  vehicles  can  be  rarefied  due  to 
very  small  length  scales,  within  a  mostly  continuum  flow  field.  For  example,  Fig.  1  shows  a  Schlieren  image  of 
typical  flow  over  a  blunted,  hypersonic  vehicle.  In  general,  the  flow  length  scales  within  the  thin  bow  shock  and 
boundary  layer  are  small  and  the  collision  length  scales  in  the  wake  are  sufficiently  large  such  that  a  microscopic 
analysis  of  these  regions  is  required  to  maintain  the  physical  accuracy  of  the  simulation.  As  the  vehicle  spends 
a  larger  portion  of  its  time  in  these  partially-rarefied  flows  with  recently  proposed  trajectory  maneuvers  [12], 
the  ability  to  maintain  physically  accurate  predictions  becomes  even  more  important  to  minimize  the  predictive 
uncertainty  and  increase  vehicle  reliability.  Recently,  a  skip  trajectory,  where  an  entry  vehicle  “skips”  through 
the  upper  atmosphere  as  illustrated  in  Fig.  3,  has  been  proposed  to  increase  the  downrange  capabilities  and 
reduce  peak  heating  of  hypersonic  entry  vehicles. 

The  AS-202  flight  case  from  the  Apollo  program  is  one  of  the  few  flight  tests  that  used  a  skip  trajectory.  Wright 
et  al.  [14]  performed  detailed  simulations  of  the  entire  flight  trajectory  using  a  continuum  description  of  the 
flow.  In  general,  they  found  that  agreement  between  available  heat  transfer  measurements  from  the  flight  tests 
and  numerical  predictions  was  good  over  most  of  the  trajectory.  However,  at  the  crest  of  the  skip  portion  of 
the  trajectory,  the  simulation  results  over-predict  the  measured  data  across  most  of  the  lee  side  of  the  aft-body. 
For  example,  Fig.  2  compares  heat  transfer  predictions  and  flight  measurements  from  a  calorimeter  located  on 
the  top  of  the  lee  side.  During  the  skip  portion,  the  CFD  simulation  method  over-predict  the  measurements  by 
over  a  factor  of  two.  Based  on  the  reconstructed  data  shown  in  Fig.  3  and  the  corresponding  Knudsen  number 
at  this  altitude  shown  in  Fig.  4,  it  is  possible  that  the  vehicle  has  transitioned  from  being  fully  continuum  to 
partially  rarefied  as  the  vehicle  has  passed  from  the  trough  to  the  crest  of  the  skip  portion  of  the  trajectory. 
Because  of  this,  Wright  et  al.  concluded  that  the  over-prediction  of  heat  transfer  may  be  due  to  the  inability 
of  the  continuum  methods  used  in  that  work  to  capture  the  important  microscopic  effects.  However,  for  this 
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Figure  1:  Schematic  of  hypersonic  flow  over  a  blunted  body  with  regions  that  typically  exhibit  non-continuum 
effects  [13] 

particular  flow,  application  of  kinetic  methods  to  the  entire  flow  would  be  computationally  expensive.  Instead, 
a  multi-scale  approach  that  only  uses  the  expensive,  microscopic  description  in  required  regions,  while  using  the 
continuum  description  throughout  the  rest  of  the  flow  field,  is  more  suitable  to  examine  any  rarefied  flow  effects. 

Simulation  of  other  multi-scale,  hypersonic  flows  is  of  particular  importance  for  developing  technologies  that 
will  enable  potential  high-mass,  Mars  missions.  Compared  to  the  Earth’s  atmosphere,  Mars’  relatively  thin 
atmosphere  causes  entry  vehicles  to  decelerate  at  much  lower  altitudes.  Depending  on  the  mass,  size,  and  shape 
of  a  vehicle,  it  may  never  reach  a  subsonic  terminal  velocity.  This  necessitates  additional  technology  to  slow 
the  vehicle  so  that  it  can  reach  subsonic  speeds  with  sufficient  time  to  prepare  for  landing.  Many  of  these  new 
technologies  require  successful  prediction  of  multi-scale  interaction  flows  to  ensure  reliability.  In  order  to  increase 
the  frontal  area  of  the  vehicle,  an  inflatable  aerodynamic  decelerator  (IAD)  (examples  shown  in  Fig.  6),  often 
called  a  ballute,  may  be  used  to  slow  the  vehicle  down  at  higher  altitudes  [15].  The  dynamic  fluid  structure 
interaction  at  deployment  will  most  certainly  occur  in  the  rarefied  regime  [16]  and  will  require  a  detailed  fluid 
structure  interaction  with  accurate  prediction  of  surface  properties  throughout  the  entire  trajectory.  In  addition, 
the  flow  field  around  the  deployed  ballute  in  the  upper  atmosphere  will  include  both  large  regions  of  continuum 
flow,  due  to  the  very  large  body,  and  small  nonequilibrium  interactions  around  tethers,  connections,  and  the 
low  density  wake  which  may  require  kinetic  analysis.  Another  possible  option  that  will  enable  high-mass,  Mars 
missions  is  supersonic  retro-propulsion  [17],  where  a  jet  is  aimed  out  the  front  of  the  aeroshell  displacing  the 
shock  and  increasing  the  total  axial  force  on  the  vehicle.  At  high  altitudes,  these  flows  may  include  a  very  dense, 
continuum  region  within  the  core  of  the  jet,  while  much  of  the  flow  may  be  rarefied.  For  example,  Fig.  7  shows 
the  flow  around  a  70°  blunted  sphere  aeroshell  with  a  propulsive  decelerator  jet  activated  at  the  center-line. 
The  collision  length  scales  within  the  jet  are  extremely  small,  such  that  the  Navier-Stokes  equations  are  valid, 
and  application  of  kinetic  analysis  to  this  region  would  be  computationally  intensive.  Conversely,  the  flow  length 
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Figure  2:  Heat  transfer  predictions  made  by  macroscopic  methods  (CFD)  on  the  after-body  of  the  AS-202  flight 
capsule  and  comparison  with  experimental  flight  data  [14] 


Figure  3:  Reconstruction  of  the  AS-202  flight  trajectory  [14] 
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Figure  4:  Variation  of  Knudsen  number  throughout  the  Earth’s  atmosphere  for  Orion  and  Apollo  capsules 


Figure  5:  Variation  of  the  gradient-length  Knudsen  number  for  flow  over  the  AS-202  capsule  at  two  different 
trajectory  points  [14] 
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Figure  6:  Typical  Inflatable  Aerodynamic  Decelerators  (IAD)  [18] 


scales  around  much  of  the  aeroshell  are  on  the  same  order  of  magnitude  as  the  collision  length  scales,  and  use  of  a 
microscopic  method  in  this  region  is  vital  to  maintain  the  necessary  physical  accuracy.  Multi-scale  methods  that 
can  accurately  simulate  these  flow  fields  will  be  an  indispensable  tool  in  development  of  these  new  technologies 
that  enable  future  scientific  and  exploration  missions  [15]. 


2.3  Plumes 

Another  class  of  multi-scale  flows  of  interest  is  rocket  exhaust  plumes.  Efficient  and  accurate  predictions  of 
atmospheric  exhaust  plumes  at  high  altitudes  are  necessary  to  ensure  that  the  chemical  rocket  maintains  efficiency 
while  also  assuring  that  the  vehicle  heating  due  to  the  complicated  flow  structure  related  to  the  under-expanded 
supersonic  jet  is  well  understood.  For  example,  Fig.  8  shows  the  continuum-rarefied  interfaces  for  flow  around  a 
sounding  rocket  traveling  in  the  upper  atmosphere.  Here,  the  flow  inside  the  nozzle  and  the  core  of  the  plume 
are  continuum,  while  the  flow  around  the  vehicle  and  in  the  outer  plume  is  rarefied.  Figure  9  shows  the  same 
rocket  during  a  simulated  stage  separation.  Now  the  flow  around  the  fore-body  of  the  separated  stage  is  entirely 
continuum,  while  the  external  flow  is  rarefied.  Application  of  full  DSMC  to  these  partially  rarefied  flow  fields 
will  require  a  large  computational  expense.  Instead,  application  of  a  multi-scale  method  will  alleviate  the  cost 
of  resolving  flow  features  in  the  highly  collisional  nozzle  and  jet  core  regions. 

In  space  flight,  knowledge  of  plume  impingement  on  surfaces  is  necessary  to  characterize  the  interaction  with 
nearby  space  craft  and  possible  deterioration  of  solar  arrays.  In  addition,  space  vehicles  landing  on  dusty  asteroids 
or  the  Moon  can  create  dust  and  debris  hazards  for  themselves  and  other  nearby  vehicles.  Accurate  prediction  of 
the  nearby  flow  fields  at  all  thrust  conditions  is  necessary  to  quantify  risk  and  develop  mitigation  strategies.  For 
example,  Lumpkin  et  al.  [20]  performed  zonally  decoupled  DSMC-CFD  simulation  of  the  RCS  plume  impingement 
of  the  reaction  control  systems  of  the  shuttle  orbiter  during  a  docking  maneuver  with  the  International  Space 
Station.  Figure  10  shows  a  planar  slice  of  the  DSMC  results  with  the  continuum  zones  highlighted.  This  work 
also  included  a  study  of  the  interaction  of  the  Apollo  Lunar  Module  jet  impinging  on  the  lunar  surface  and  the 
formation  of  a  crater  of  dust.  An  illustration  of  a  typical  simulation  is  shown  in  Fig.  11.  Their  study  results 
were  mixed  and  found  that  the  decoupled  approach  may  increase  the  computational  cost  compared  to  a  coupled 
approach  since  the  CFD-DSMC  interface  region  had  to  be  placed  well  into  the  continuum  regime  to  be  upstream 
of  any  impingement  shock  waves.  This  resulted  in  DSMC  only  resolving  gross  flow  features,  while  finer  features, 


RTO-EN-AVT-194 


6-7 


Hybrid  Particle-Continuum  Numerical  Methods  for  Aerospace  Applications 


ORGANIZATION 


-0.005  - 


-0.01  - 


-0.01  -0.005  0 


0.005 

X  [m] 


o.oi 


0.015 


0.02 


Figure  7:  Continuum/rarefied  domains  for  flow  over  an  aeroshell  with  a  propulsive  decelerator  activated  in  a 
partially  rarefied  flow  field  with  a  nominal  global  Knudsen  number  of  0.002  based  on  the  aeroshell  diameter  and 
free  stream  Mach  number  of  12 


Figure  8:  Continuum/rarefied  domains  for  a  high  altitude  rocket  plume  [19] 
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Figure  9:  Continuum/rarefied  domains  for  a  rocket  at  stage  separation  [19] 


such  as  surface  shear,  have  visible,  unphysical  numerical  artifacts.  Instead,  a  fully  coupled  hybrid  method  may 
provide  a  significant  reduction  in  computational  expense,  with  little  loss  of  physical  accuracy. 


2.4  Viscous  Interaction  Flows 

Hypersonic  viscous  interaction  flows  involve  the  interactions  between  shock  waves  and  boundary  layers,  producing 
complex  and  highly  non-linear  flow  fields.  Shock-boundary  layer  interactions  occur  in  inlets  of  supersonic  and 
hypersonic  propulsion  devices  and  in  the  vicinity  of  control  surfaces  on  hypersonic  vehicles.  In  such  flows,  shock- 
shock  interactions  often  produce  reflected  shock  waves  that  impinge  on  the  surface  of  the  vehicle.  The  interaction 
between  a  strong  shock  wave  and  a  boundary  layer  often  causes  the  flow  to  separate  and  form  a  region  of  high- 
pressure,  recirculating  gas  next  to  the  surface.  High-speed  flow  hitting  such  a  separation  region  significantly  alters 
shock  structures  and  interaction  regions  which,  in  turn,  affect  the  extent  of  flow  separation.  Peak  aerothermal 
loads  are  observed  at  the  location  of  shock  impingement  on  the  vehicle  surface  and  therefore  accurate  prediction 
of  this  phenomenon  is  important  in  the  design  of  a  hypersonic  vehicle.  In  addition,  the  heating  and  frictional 
loads  generated  within  a  separated  region  are  drastically  different  than  when  the  flow  is  purely  attached.  Such 
differences  may  reduce  the  effectiveness  of  a  control  surface  and  thus  accurately  predicting  the  extent  of  any  flow 
separation  is  also  a  very  important  aspect  in  the  design  of  hypersonic  vehicles. 

An  example  of  a  viscous  interaction  flow  is  shown  in  Fig.  12  which  was  part  of  a  blind  validation  study  with 
experimental  measurements  taken  in  the  Large  Energy  National  Shock  (LENS)  facility  at  the  Calspan  -  University 
of  Buffalo  Research  Center  (CUBRC)  [21].  For  Run  11,  which  had  a  global  Knudsen  number  of  0.0008,  it  has 
been  found  that  solutions  to  the  Navier-Stokes  equations  predict  a  separation  bubble  that  forms  along  the 
surface  that  is  significantly  larger  than  experimental  measurements.  In  general,  DSMC  predictions  of  the  same 
case  have  better  agreement  with  the  experimental  measurements  in  this  leading  edge  region,  but  require  a  large 
computational  expense  to  accurately  resolve  the  shock  impingement  downstream  of  the  flare  junction  [21,  22]. 
Instead,  if  a  multi-scale  method  can  apply  a  continuum  solver  in  regions  where  the  mean  free  path  is  small,  and 
the  Navier-Stokes  equations  are  valid,  the  computational  requirements  of  obtaining  a  physically  accurate  solution 
may  be  reduced. 
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Figure  10:  Plume  impingement  study  of  the  shuttle  orbiter  docking  with  the  international  space  station  [20] 
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Figure  11:  Plume  impingement  study  of  the  Apollo  Lunar  Module  with  the  lunar  surface  [20] 
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Figure  12:  Geometry  and  variation  in  mean  free  path  for  a  hypersonic  shock-boundary  layer  interaction  test  case 
[23] 

3  Hybrid  Particle-Continuum  Framework 

This  section  provides  an  overview  of  typical  challenges  associated  with  the  development  of  hybrid  particle- 
continuum  simulation  techniques  along  with  published  solutions.  First,  an  overview  of  methods  to  compute  the 
demarcation  between  continuum  and  rarefied  regions  of  the  flow  is  provided.  Then,  typical  coupling  procedures 
used  to  transfer  information  between  flow  modules  is  outlined. 

3.1  Continuum  Breakdown 

Both  the  accuracy  and  efficiency  of  a  hybrid  DSMC-CFD  method  depend  strongly  on  proper  placement  of  the 
interface  location.  For  physical  accuracy,  the  interface  must  be  located  within  regions  that  can  be  considered 
in  collisional  equilibrium,  where  the  velocity  distribution  is  only  slightly  perturbed  from  equilibrium,  and  the 
Navier-Stokes  equations  are  valid.  To  maximize  efficiency,  a  hybrid  method  requires  the  interface  between  CFD 
and  DSMC  to  be  located  near  the  edge  of  the  collisional  equilibrium  region.  Empirical  breakdown  parameters 
have  been  used  both  as  switching  criteria  for  hybrid  rarefied-continuum  methods  and  also  in  analysis  of  continuum 
predictions  to  ensure  that  the  Navier-Stokes  equations  are  valid  throughout  the  flow  field. 

Bird  [24]  studied  DSMC  simulations  of  expansion  flows  and  proposed  a  continuum  breakdown  parameter,  P,  for 
expansions  shown  in  Eq.  2  where  v  is  the  collision  frequency  and  p  is  the  mass  density.  Through  comparison 
of  many  one  dimensional  simulations  of  spherical  and  cylindrical  expansions  with  the  DSMC  method,  Bird 
found  that  breakdown  of  translational  equilibrium  occurred  when  P  exceeded  0.05.  Bird  also  showed  that  this 
parameter  can  be  used  to  estimate  the  breakdown  of  rotational  equilibrium.  However,  the  simulation  methods 
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used  a  constant  rotational  collision  number  of  2,  which  under  predicts  the  rotational  relaxation  process  for  typical 
gases  [25]. 

v 


D  (In  p) 


Dt 


(2) 


Boyd  et  al.  [26,  27]  have  proposed  a  gradient-length  Knudsen  number,  shown  in  Eq.  3  where  A  is  the  local  mean 
free  path  and  Q  is  some  macroscopic  flow  quantity  of  interest,  such  as  density,  flow  speed,  or  temperature. 
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The  local  mean  free  path  used  in  the  gradient-length  Knudsen  number  is  calculated  using  Eq.  4  where  n  is  the 
number  density,  Ttra  is  the  translational  temperature,  u>  is  the  macroscopic  viscosity  temperature  exponent, 
and  Tref  is  the  temperature  that  the  reference  cross  section,  oye/,  is  calculated  at,  which  is  consistent  with 
the  variable  hard  sphere  collision  model  [1]  used  in  the  DSMC  simulation  and  the  corresponding  temperature 
viscosity  relation  used  within  the  CFD  module.  The  reference  cross  section  can  be  calculated  using  Eq.  5  where 
k  is  the  Boltzmann  constant  and  pref  is  the  coefficient  of  viscosity  at  Tref. 


Tref 
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Gref  ~  2  (5  -  2 uj)  (7  -  2uj)  pref  ^5) 

The  gradient-length  Knudsen  numbers  are  directly  related  to  the  Chapman-Enskog  expansion  terms  for  mass, 
momentum,  and  energy  molecular  diffusion.  Through  detailed  comparison  of  kinetic  and  continuum  simulations 
of  planar  shock  waves,  spherically  blunted  bodies,  and  interaction  flows,  Boyd  et  al.  found  that  in  regions 
where  the  maximum  gradient-length  Knudsen  number  exceeded  0.05,  the  difference  between  full  DSMC  and 
full  CFD  exceeds  5%.  Therefore,  a  conservative  breakdown  parameter  of  0.05,  such  that  DSMC  is  used  when 
max  (-KtzGIj_q)  >  0.05  and  CFD  is  used  elsewhere,  should  not  introduce  more  than  5%  error  within  a  hybrid 
method.  Recent  results  [28,  29]  with  a  hybrid  DSMC-CFD  method  have  shown  that  relaxing  the  breakdown 
parameter  to  a  value  of  0.1  still  produces  hybrid  predictions  that  remain  in  excellent  agreement  with  full  DSMC 
simulation  results. 

Similarly,  Garcia  et  al.  [30]  have  proposed  a  breakdown  parameter  based  directly  on  the  expansion  terms  of  the 
Chapman-Enskog  distribution  function,  which  is  shown  in  Eq.  6  and  where  the  normalized  shear  stress  and  heat 
flux  terms  are  shown  in  Eqs.  7  and  8,  respectively.  A  conservative  switching  criterion  of  0.1,  such  that  DSMC  is 
used  when  B  >  0.1  and  CFD  is  used  otherwise,  has  been  proposed  [31]  while  in  practice,  the  switching  criteria 
has  been  successfully  relaxed  to  0.2  for  certain  flows  [32]. 


B  =  max  ,  |g*|) 

(6) 

/i  /  dui  duj  2  duk  -  A 
p  \cAcj  ^  dxi  3  dxk  ’  / 

(7) 

*  k  j~2m  dT 

(8) 

^  p\kT  dxi 

6  - 12 


RTO-EN-AVT-194 


Hybrid  Particle-Continuum  Numerical  Methods  for  Aerospace  Applications 


Tiwari  [33]  has  proposed  a  switching  criteria,  ||<^||,  as  shown  in  Eq.  9  where  q  is  the  heat  flux  vector,  r  is  the 
shear  stress  tensor,  p  is  the  local  mass  density,  R  is  the  universal  gas  constant,  and  T  is  the  local  translational 
temperature.  For  the  continuum  approximation  to  hold,  Tiwari  reasoned  that  this  parameter  must  be  much  less 
than  unity,  therefore  Tiwari  stated  that  kinetic  solvers  should  be  used  anywhere  that  ||</>||  >  e,  where  e  is  much 
less  than  unity. 


2m! 

5  RT 
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Recently,  Lockerby  et  al.  [34]  have  proposed  a  switching  criterion,  shown  in  Eq.  10,  that  compares  the  difference 
between  typical  shear  stress  and  heat  flux  terms  used  to  close  the  Navier-Stokes  equations  (Newtonian  fluid  and 
Fourier’s  Law)  and  the  third  order  physically  accurate  Regularized  13-Moment  (R13)  equations  [35].  They  have 
also  formulated  an  estimate  of  the  R13  shear  stress  and  heat  flux  terms  using  only  information  available  from 
the  Navier-Stokes  equations  so  that  the  breakdown  parameter  can  be  calculated  without  a  full  R13  solution. 
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In  addition  to  collisional  nonequilibrium,  a  hybrid  method  must  ensure  that  DSMC  is  used  in  other  regions  of 
the  flow  where  certain  physical  processes,  such  as  rotational  relaxation,  are  important  but  are  ignored  in  the 
continuum  module.  For  example,  if  the  continuum  module  assumes  that  translational  and  rotational  modes  are  in 
equilibrium,  then  the  hybrid  method  should  ensure  that  regions  that  can  be  considered  in  collisional  equilibrium, 
such  that  the  velocity  distribution  function  can  be  described  by  a  Chapman-Enskog  velocity  distribution  function, 
but  where  the  rotational  mode  is  not  in  equilibrium  with  the  translational  mode  be  assigned  to  the  DSMC  module. 
For  near  equilibrium  flows  over  blunt  bodies,  which  are  of  interest  for  hybrid  methods,  almost  the  entire  region 
behind  the  strong  expansion  displays  significant  nonequilibrium  between  the  rotational  and  translational  energy 
modes.  This  would  greatly  increase  the  size  of  the  region  simulated  with  the  DSMC  module  in  a  hybrid  simulation 
and  have  a  serious  adverse  effect  on  the  numerical  efficiency.  Schwartzentruber  et  al.  [36,  28,  37]  found  that  only 
the  strong  thermal  relaxation  behind  a  strong  bow  shock  had  a  large  adverse  effect  on  the  physical  accuracy  of 
the  hybrid  method  and  proposed  an  additional  breakdown  parameter  shown  in  Eq.  11  with  a  switching  parameter 
of  0.01  so  that  regions  where  the  translational  temperature  exceeds  the  rotational  temperature  by  more  than  1% 
are  simulated  with  the  DSMC  method. 
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In  addition  to  the  strong  bow  shock,  this  additional  breakdown  parameter  has  been  found  to  increase  the  size  of 
the  DSMC  region  behind  the  recompression  wave  in  the  wake  of  blunt  body  flows  [38].  This  has  an  adverse  effect 
on  the  efficiency  of  a  hybrid  method,  but  is  required  to  maintain  sufficient  physical  accuracy.  By  inclusion  of  a 
separate  rotational  energy  equation  in  the  continuum  module,  nearly  the  entire  region  that  displays  rotational 
nonequilibrium  with  near  translational  equilibrium  can  be  instead  simulated  within  the  continuum  module  with 
little  physical  error. 


However,  even  if  a  separate  rotational  energy  equation  is  included  in  the  continuum  module,  the  additional 
breakdown  parameter  is  still  needed  to  ensure  the  the  rotational  energy  probability  density  function  is  near 
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equilibrium  in  the  continuum  region,  but  the  breakdown  value  can  be  relaxed  to  0.1  such  that  DSMC  is  used  in 
all  regions  where  the  rotational  temperature  deviates  from  the  translational  temperature  by  more  than  10%. 

For  example,  Fig.  13  shows  a  comparison  of  translational  temperature  prediction  of  a  hybrid  method  with  (top) 
and  without  (bottom)  the  rotational  nonequilibrium  switch  activated  along  with  full  DSMC  simulation  results. 
Without  the  nonequilibrium  switch,  the  interface  location  between  the  CFD  and  DSMC  modules  is  located 
very  near  the  shock.  Since  CFD  does  not  contain  the  physical  accuracy  required  to  model  this  portion  of  the 
rotational  relaxation  process  as  the  rotational  energy  distribution  function  is  highly  non-Boltzmann  in  this  region, 
the  post  shock  temperatures  are  over  predicted  compared  to  full  DSMC.  With  a  rotational  nonequilibrium  cutoff 
parameter  of  0.1,  such  that  DSMC  is  used  in  regions  where  the  difference  between  the  translational  and  rotational 
temperatures  exceeds  10%  of  the  rotational  temperature,  the  portion  of  the  flow  field  modeled  with  the  continuum 
solver  is  slightly  decreased  to  only  regions  where  the  continuum  solver  is  physically  valid  and  agreement  between 
the  hybrid  method  and  full  DSMC  results  is  greatly  improved. 

In  addition,  velocity  and  rotational  energy  distribution  functions  are  sampled  from  the  full  DSMC  solution  at 
locations  at  the  edge  of  the  interface  location  predicted  by  each  of  the  hybrid  simulations  denoted  as  A  and 
B  in  Fig.  13.  Figures  14a  and  14b,  respectively,  compare  the  velocity  distribution  functions  and  rotational 
energy  distribution  function  sampled  at  point  A  from  DSMC  with  equilibrium  theory.  Although  the  velocity 
probability  density  functions  are  in  excellent  agreement  with  equilibrium  theory,  the  rotational  energy  distribution 
function  differs  significantly  from  equilibrium.  The  peak  in  the  rotational  energy  probability  density  function 
at  low  rotational  energies  signifies  that  there  is  still  a  large  number  of  particles  that  have  not  experienced  a 
rotationally  inelastic  collision  at  the  post  shock  temperature  and  still  maintain  the  slope  of  an  equilibrium 
probability  density  function  at  the  free  stream  temperature.  In  contrast,  Figs.  15a  and  15b  compare  the  velocity 
and  rotational  energy  distribution  functions,  respectively,  at  the  continuum  interface  location  computed  with 
the  added  breakdown  parameter  at  point  B.  At  this  point  in  the  flow,  both  the  velocity  and  rotational  energy 
probability  density  functions  are  in  much  better  agreement  with  the  equilibrium  description  and  the  models 
used  in  the  continuum  solver  are  valid.  Since  the  equilibrium  rotational  energy  distribution  function  is  calculated 
based  on  the  average  rotational  energy,  comparison  of  higher  order  moments,  such  as  the  variance,  calculated  from 
sampled  data  and  the  equilibrium  distribution  can  be  used  as  a  measure  of  degree  of  rotational  nonequilibrium. 
At  point  A,  the  sampled  variance  differs  by  nearly  25%  from  equilibrium,  while  the  sampled  variance  differs  by 
less  than  10%  from  the  equilibrium  value  at  point  B. 

3.2  Coupling  Procedures 

The  method  of  information  transfer  between  particle  and  continuum  methods  can  typically  be  split  into  two 
categories:  coupling  by  maintaining  consistent  fluxes  or  coupling  by  maintaining  consistent  state  properties 
in  reservoir  cells.  Figure  16  provides  a  visual  comparison  of  the  two  coupling  methods.  A  flux-based,  coupling 
method,  which  is  depicted  in  Fig.  16(a),  involves  calculating  fluxes  of  conserved  variables  at  the  interface  location 
according  to  particle  cell,  FP,  and  according  to  the  continuum  cell,  Fc.  The  particle  flux  can  be  directly  calculated 
by  tracking  the  number  of  particles  that  cross  the  interface,  but  the  continuum  flux  must  be  extrapolated  using 
macroscopic  cell-averaged  state  quantities  and  their  gradients.  In  general,  the  estimated  fluxes  are  not  the 
same  and  are  often  modified  to  ensure  that  mass,  momentum,  and  energy  are  conserved  within  the  simulation. 
This  modified  flux  is  applied  as  a  boundary  condition  to  both  simulation  techniques  and  is  used  to  update  the 
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Figure  13:  Comparison  of  translational  temperatures  predicted  by  the  MPC  method  with  (top)  and  without 
(bottom)  the  rotational  nonequilibrium  breakdown  parameter  compared  with  full  DSMC  of  Mach  12  flow  over  a 
cylinder  with  a  global  Knudsen  number  of  0.01  [38] 
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Figure  14:  Comparison  of  probability  density  functions  predicted  by  DSMC  and  equilibrium  theory  at  point  A 
shown  in  Figure  13  [38] 
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Figure  15:  Comparison  of  probability  density  functions  predicted  by  DSMC  and  equilibrium  theory  at  point  B 
shown  in  Figure  13  [38] 
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(a)  Flux-based  coupling.  (b)  State-based  coupling. 


Figure  16:  Schematic  of  typical  hybrid  coupling  procedures 

continuum  solution  and  to  estimate  the  probability  density  functions  that  are  sampled  to  assign  DSMC  particle 
information. 

With  state-based  coupling,  shown  in  Fig.  16(b),  transfer  of  information  between  simulation  methods  is  performed 
by  one  method  providing  a  Dirichlet  boundary  condition  to  the  other  simulation  technique.  Each  method  can 
then  update  its  respective  regions  using  algorithms  internal  to  the  corresponding  simulation  method.  For  coupling 
from  the  particle  to  the  continuum  method,  these  values  are  calculated  by  performing  an  averaging  procedure  to 
sample  state  quantities  of  interest,  such  as  average  density,  velocity,  and  energy  in  DSMC  cells  along  the  edge 
of  the  interface  location.  These  values  are  assigned  to  continuum  ghost  cells  which  are  used  to  calculate  inviscid 
and  viscous  fluxes  to  update  the  solution  within  the  continuum  domain.  For  coupling  from  the  continuum  to 
the  particle  method,  average  state  and  gradient  information  is  used  to  estimate  the  probability  density  functions 
in  DSMC  boundary  cells.  Before  each  particle  iteration,  all  simulators  in  these  boundary  cells  are  deleted  and 
new  particles  are  generated  consistent  with  the  estimated  probability  density  functions  constructed  from  the 
continuum  data.  Then,  simulators  in  reservoir  cells  are  allowed  to  move  throughout  the  domain  using  the 
internal  DSMC  move  algorithm. 

For  both  coupling  procedures,  techniques  to  reduce  the  statistical  scatter  are  often  employed  to  ensure  that 
the  continuum  method  remains  stable.  For  many  flows  that  involve  strong  two-way  coupling,  the  efficiency  of 
the  hybrid  method  is  determined  by  how  often  information  is  passed  between  flow  modules.  Unfortunately,  the 
statistical  scatter  of  averaged  DSMC  data  over  very  few  time-steps  can  have  large  variations  at  the  continuum- 
rarefied  boundary  that  could  create  numerical  instabilities  in  the  continuum  module.  Many  techniques  to  reduce 
this  statistical  scatter  have  been  proposed  for  different  hybrid  methods  and  will  be  reviewed  in  the  following 
subsection. 

Often,  an  overlap  region  between  the  continuum  and  particle  boundary  cells  is  also  used,  where  the  particle 
region  is  extended  into  the  continuum  domain  and  both  simulation  methods  calculate  the  solution.  These  overlap 
cells  can  be  used  to  filter  the  physical  inaccuracies  from  an  incorrect  initial  solution  or  to  compare  simulation 
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predictions  to  ensure  that  the  interface  location  is  located  in  a  near  continuum  region,  where  the  Navier-Stokes 
equations  are  appropriate. 


4  Hybrid  Codes  and  Their  Application 

Over  the  past  two  decades,  many  studies  of  hybrid  particle-continuum  techniques  have  been  performed  that  have 
continuously  expanded  the  state  of  the  art  capabilities. 

Wadsworth  and  Erwin  developed  a  strongly  coupled,  flux-based  hybrid  DSMC-CFD  scheme  and  applied  it  both 
to  simulate  one  dimensional  shocks  [39]  and  two  dimensional  rarefied  slit  flow  [40].  In  these  studies,  a  Maxwellian 
distribution  was  used  to  generate  simulation  particles  at  the  interface  location,  which  was  consistent  with  the 
governing  Euler  equations  in  the  continuum  region  while  the  domain  boundaries  remained  fixed.  Although  the 
flow  was  unsteady,  the  scheme  used  cumulative  averages  to  calculate  the  particle  to  continuum  flux  to  reduce 
the  inherent  statistical  scatter,  while  the  continuum  to  particle  flux  was  time  accurate  based  on  the  current 
continuum  solution. 

Hash  and  Hassan  also  coupled  particle  and  continuum  codes  using  a  flux-based  procedure  to  simulate  Couette 
flow  [41]  and  near  equilibrium,  hypersonic  flow  over  a  sting-mounted  blunted  sphere-cone  [42,  43].  Bird’s  DSMC 
method  that  took  into  account  internal  excitation  and  finite  rate  chemistry  was  used  as  the  DSMC  module 
while  the  modified  Navier-Stokes  equations  were  solved  on  a  structured  grid  in  the  continuum  module.  Both  the 
Marshak  condition  and  property  interpolation  technique  were  employed  to  transfer  information  at  the  interface. 
Along  much  of  the  rarefied-continuum  interface  location  for  the  hypersonic  flow  case,  the  normal  Mach  number 
was  small  which  greatly  increased  the  statistical  scatter  associated  with  particle  fluxes  at  the  boundary.  The 
statistical  scatter  was  reduced  using  a  smoothing  operator  over  highly  sampled  data  for  the  property  extrapolation 
before  imposing  the  flux  boundary  condition.  Although  these  smoothing  procedures  reduced  the  scatter  for 
bulk  properties,  oscillations  of  properties  with  large  gradients  along  the  boundary  could  produce  unphysical, 
negative  values  of  density.  In  addition,  the  statistical  scatter  of  gradients  were  still  appreciable  with  low  order 
smoothing  operators  so  a  fourth-order  smoothing  procedure  with  a  7-point  stencil  had  to  be  used  in  certain 
regions.  To  maintain  consistency  at  the  interface,  the  researchers  also  found  that  sampling  particle  velocities 
from  the  Chapman-Enskog  velocity  probability  density  function  was  necessary  in  regions  that  displayed  significant 
momentum  or  energy  transfer  [41]. 

Roveda  et  al.  proposed  a  hybrid  particle-continuum  method  that  used  a  state-based  coupling  procedure  to 
simulate  moving  planar  shock  waves  [44]  and  two  dimensional  unsteady  flow  [45].  The  continuum  module  solved 
the  Euler  equations  using  Nadiga’s  adaptive  discrete  velocity  (ADV)  method  [46],  while  Bird’s  DSMC  method 
was  used  in  the  particle  module.  Although  the  state-based  procedure  had  less  statistical  scatter  at  the  boundary 
when  compared  to  flux-based  coupling  procedures,  the  method  was  time-accurate  which  reduced  the  number 
of  statistically  independent  samples  to  the  current  number  of  DSMC  simulator  particles  in  each  cell.  In  order 
to  reduce  the  statistical  scatter  in  boundary  cells,  their  method  used  a  novel  algorithm  to  clone  particles  in  a 
buffer  around  the  rarefied  region  to  increase  the  number  of  nearly  statistically  independent  DSMC  samples.  This 
had  an  effect  of  reducing  the  statistical  scatter  by  over  a  factor  of  two.  Two  layers  of  DSMC  reservoir  cells 
were  constructed  along  the  DSMC  boundary  and  particles  were  assigned  velocities  from  the  Maxwell-Boltzmann 
velocity  distribution  function  which  is  consistent  with  the  continuum  solver  used  in  this  work. 
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Garcia  et  al.  proposed  an  adaptive  mesh  and  algorithm  refinement  method  [31,  32]  that  introduces  the  DSMC 
methodology  at  the  highest  level  of  refinement  of  an  adaptive  mesh  scheme.  During  each  time-step,  the  continuum 
method  is  applied  to  the  entire  flow.  Then,  regions  that  are  labeled  as  nonequilibrium  regions  are  resolved  with 
a  finer  mesh  that  meets  DSMC  requirements  and  multiple  DSMC  time-steps  are  taken  to  match  the  larger, 
continuum  time-step.  Particles  are  generated  in  reservoir  cells  around  DSMC  regions  consistent  with  the  time- 
interpolated  continuum  data,  while  particle  fluxes  across  each  continuum  face  are  recorded.  After  the  DSMC 
simulation  procedure  has  reached  the  current  continuum  time,  DSMC  samples  are  used  to  update  the  applicable 
continuum  cells  and  conserved  fluxes  are  calculated  based  on  the  sampled  particle  simulator  fluxes.  This  method 
has  been  successfully  applied  to  one-  and  two-dimensional  gas  flows. 

Sun  et  al.  [47]  have  proposed  a  hybrid  method  that  couples  a  Navier-Stokes  solver  with  a  novel  information 
preservation  (DSMC-IP)  method  to  simulate  micro-scale  flows.  The  IP  method  is  an  extension  of  Bird’s  DSMC 
method  that  tracks  additional  macroscopic  variables  in  each  DSMC  cell  that  contain  less  statistical  scatter.  The 
reduced  statistical  scatter  of  the  IP  method  allows  it  to  be  tightly  coupled  with  a  Navier-Stokes  solver  with  little 
effect  on  the  stability  in  continuum  regions.  They  found  that  Garcia’s  breakdown  parameter  [30,  31]  worked  best 
for  near  equilibrium,  micro-scale  gas  flows.  Figure  17  shows  a  schematic  of  the  state-based  coupling  procedure  at 
the  continuum-rarefied  interface  location.  Reservoir  cells  are  used  and  particle  velocities  are  assigned  by  sampling 
the  Chapman-Enskog  velocity  distribution  function  using  the  continuum  information  in  the  corresponding  cells. 
Application  of  the  NS/DSMC-IP  hybrid  method  has  been  successful  on  many  micro-scale  flows  including  Couette 
flow,  Rayleigh  flow,  and  external  flow  over  micro-airfoils.  For  example,  Fig.  18  shows  a  comparison  of  flow  velocity 
around  a  flat  plate  at  zero  angle  of  attack  predicted  by  full  Navier-Stokes,  full  DSMC-IP,  and  hybrid  NS/DSMC- 
IP.  The  hybrid  method  has  been  able  to  successfully  reproduce  microscopic  effects  in  the  rarefied,  near-body 
region  that  the  continuum,  Navier-Stokes  equations  are  unable  to  capture,  and  has  been  used  to  characterize 
external  flow  over  micro-plates  at  various  angles  of  attack  and  Reynolds  numbers.  For  example,  Fig.  19  shows 
a  hybrid  NS/DSMC-IP  prediction  of  flow  over  a  flat  plate  with  5%  thickness  at  a  20°  angle  of  attack.  Here  the 
rarefied  region,  which  is  highlighted  in  Fig.  20,  is  confined  to  a  small  region  near  the  leading  edge.  Rather  than 
simulate  the  entire  flow  with  a  rarefied  flow  method,  this  hybrid  method  is  able  to  apply  the  DSMC-IP  method 
to  only  regions  that  exhibit  translational  nonequilibrium  effects. 

Wang  and  Boyd  [48]  have  proposed  using  the  DSMC-IP  method  coupled  to  a  Navier-Stokes  solver  within  a  hybrid 
framework  to  simulate  partially  rarefied,  hypersonic  flow.  Although  the  method  was  successful  for  certain  two- 
dimensional  flows,  it  was  found  that  the  DSMC-IP  scheme  produced  an  incorrect  post-shock  state  and  shock  wave 
thickness  when  applied  to  planar  shock  waves.  A  new  formulation  for  the  DSMC-IP  energy  flux  [49]  provided 
a  partial  remedy  to  the  discrepancies,  but  at  a  large  computational  expense.  The  hybrid  NS/DSMC-IP  method 
has  been  applied  to  planar  shock  waves,  a  blunted  cone  with  a  small  nose  radius,  and  the  hollow  cylinder-flare 
experimental  case  that  was  described  in  the  viscous  interaction  flow  subsection. 

Wu  et  al.  [50]  have  proposed  a  parallel  coupled  DSMC-NS  method  to  simulate  hypersonic  flow  over  a  25°  wedge 
and  expansion  of  nitrogen  gas  from  a  three-dimensional  nozzle  into  a  near-vacuum  container.  The  authors 
used  a  shell  script  that  coupled  existing  DSMC  and  CFD  codes.  State-based  coupling  was  used  to  transfer 
information  between  the  two  simulation  methods.  The  DSMC  region  is  extended  into  the  continuum  region  and 
the  breakdown  parameter  is  continuously  applied  to  the  flow  field  to  ensure  that  the  interface  location  remains 
in  a  continuum  location.  Though  the  continuum  method  solved  the  Navier-Stokes  equations,  the  particles 
generated  at  the  continuum-rarefied  interface  were  assigned  velocities  consistent  with  the  Maxwell-Boltzmann 
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Figure  17:  Schematic  of  DSMC-IP/NS  coupling  [47] 
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Figure  18:  Prediction  of  micro-scale  flow  over  a  flat  plate  with  DSMC-IP,  Navier-Stokes  solver,  and  hybrid 
DSMC-IP/NS  [47] 
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Figure  19:  Prediction  of  Mach  0.2  flow  over  a  flat  plate  at  a  20°  angle  of  attack  with  hybrid  DSMC-IP/NS  [47] 
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Figure  20:  Enlargement  of  the  continuum/rarefied  demarcation  near  the  leading  edge  of  Mach  0.2  flow  over  a 
flat  plate  at  a  20°  angle  of  attack  [47] 
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velocity  distribution  function  to  increase  efficiency.  Initially,  the  hybrid  method  used  the  maximum  gradient 
length  Knudsen  number  with  a  cutoff  parameter  of  0.02  to  ensure  that  the  interface  location  was  in  a  region  that 
was  very  near  equilibrium.  However,  when  compared  to  a  full  DSMC  simulation,  the  hybrid  simulation  took 
about  20%  more  computational  time.  The  continuum  code  used  a  single  temperature  to  describe  the  translational 
and  rotational  modes  which  required  an  additional  thermal  breakdown  parameter  very  similar  to  that  shown  in 
Eq.  11  to  ensure  that  regions  that  displayed  rotational  nonequilibrium  effects  are  simulated  with  the  DSMC 
method.  To  improve  the  efficiency,  later  work  [51]  used  only  the  gradient  length  Knudsen  number  based  on 
pressure  which  then  includes  most  of  the  boundary  layer  to  be  simulated  with  the  continuum  method  instead 
of  the  DSMC  module.  The  cutoff  parameter  was  also  relaxed  to  a  value  of  0.05.  In  addition,  the  rotational 
breakdown  parameter  was  reformulated  to  Eq.  12  where  Tq  is  the  temperature  associated  with  the  gth  energy 
mode,  is  the  number  of  rotational  degrees  of  freedom,  and  (v  is  the  number  of  vibrational  degrees  of  freedom. 
A  cutoff  parameter  value  of  0.03  was  employed  for  this  additional  switch  such  that  DSMC  is  used  in  regions 
where  P^he  >  0.03. 


D*  _ 

rThe  ~ 


\ 


3  +  Cr  +  C« 


(12) 


The  number  of  coupling  iterations  was  also  reduced  from  six  to  four  which  resulted  in  the  hybrid  method  requiring 
nearly  a  factor  of  2  less  computational  time  compared  to  the  corresponding  full  DSMC  simulation. 


Schwartzentruber  and  Boyd  [52]  have  proposed  a  Modular  Particle-Continuum  (MPC)  method  to  simulate  par¬ 
tially  rarefied,  steady-state  flow  over  hypersonic  vehicles.  This  method  couples  existing  DSMC  and  CFD  modules 
that  remain  nearly  unmodified.  State-based  coupling  is  performed  with  particle  velocities  assigned  by  sampling 
of  the  Chapman-Enskog  velocity  distribution  function  using  the  algorithm  of  Garcia  and  Alder  [30],  while  the 
statistical  scatter  of  particle  samples  is  reduced  by  use  of  a  novel  subrelaxation  scheme  of  Sun  and  Boyd  [53] 
which  is  shown  in  Eq.  13  where  Q  is  a  quantity  of  interest,  j  is  the  current  iteration,  and  <I>  is  the  subrelaxation 
parameter.  Typically  a  subrelaxation  parameter  value  of  0.001  is  used. 


(Q)j  =  (i  -  *)  (Q)j- 1  +  *Qs  (13) 

The  method  has  been  successfully  applied  to  reproduce  full  DSMC  predictions  for  planar  shock  waves  [52], 
two-dimensional  blunt  body  flows  [54,  36]  and  axi-symmetric  blunt  body  [28]  and  viscous  interaction  flows  [23]. 
Though  the  overall  hybrid  cycle  is  similar  to  that  proposed  by  Wu  et  ah,  there  is  one  subtle  and  important 
difference.  The  MPC  method  does  not  call  the  continuum  solver  until  the  interface  has  stopped  moving  to  ensure 
that  the  continuum  method  is  applied  only  in  regions  where  the  Navier-Stokes  equations  are  valid.  This  is  vital 
for  the  proper  simulation  of  the  shock  layers  over  hypersonic,  blunt  bodies  to  correctly  predict  the  post-shock 
condition  and  shock  standoff  distance  [55]. 

Recently,  extensions  of  the  MPC  method  have  allowed  for  simulation  of  rotational  [38]  and  vibrational  nonequi¬ 
librium  [56]  throughout  the  entire  flow  domain,  and  the  method  has  been  parallelized  [57]  for  distributed  memory 
systems  using  dynamic  domain  decomposition  routines.  The  computational  speed-up  of  the  MPC  method  over 
full  DSMC  has  ranged  from  about  1.5  for  flows  that  are  nearly  entirely  rarefied  to  over  a  factor  of  30  for  near 
equilibrium  flows. 
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Although  the  sub  relaxation  technique  reduces  the  scatter  of  macroscopic  flow  variables,  it  should  only  be  used 
if  the  flow  variables  of  interest  are  resolved  in  less  than  ^  iterations.  For  almost  all  flow  variables,  this  can  often 
be  done  in  very  few  iterations.  For  vibrational  energies,  a  discrete  probability  density  function  is  used  that  is 
consistent  with  the  simple  harmonic  oscillator  assumption  and  the  number  of  iterations  required  to  resolve  a 
low  vibrational  temperature  becomes  enormous.  Figure  21  compares  the  subrelaxation  average  rotational  and 
vibrational  temperatures  as  a  function  of  iteration  number  with  20  particles  per  iteration  at  various  values  of 
temperature.  In  general,  the  rotational  temperature  (and  also  translational  temperature  that  is  not  shown)  can 
be  resolved  in  very  few  iterations  and  a  subrelaxation  parameter  of  0.001  is  sufficient  regardless  of  the  mean  energy 
content  since  the  discrete  energy  steps  of  these  modes  are  much  less  than  the  mean  particle  energy  content  and 
the  energy  probability  density  functions  can  be  considered  continuous.  In  contrast,  most  free  stream  vibrational 
energies  are  much  less  than  the  discrete  energy  step  size  and  the  probability  of  experiencing  a  vibrationally  excited 
molecule  is  extremely  low.  For  example,  when  the  vibrational  temperature  is  less  than  O.10VIB  (which  corresponds 
to  a  free  stream  temperature  of  less  than  339.5  K  for  IV2),  less  than  two  vibrationally  excited  molecules  are  found 
in  each  set  of  A  iterations.  Decreasing  the  subrelaxation  parameter,  d>,  could  decrease  the  statistical  scatter 
of  the  subrelaxation  averaged  vibrational  temperature  at  the  expense  of  efficiency,  but  typical  high  altitude 
temperatures  would  require  <I>  to  be  many  orders  of  magnitude  smaller  than  what  is  currently  used.  This  would 
make  any  coupled  hybrid  method  using  this  technique  slower  than  full  DSMC. 


Figure  21:  Comparison  of  the  level  of  statistical  scatter  of  subrelaxation  averages  of  internal  temperatures  at 
various  levels. 


Instead  of  assigning  vibrational  energies  consistent  with  the  discrete  Boltzmann  energy  probability  density  func¬ 
tion,  the  average  vibrational  energy  is  assigned  to  all  particles  in  the  boundary  cells.  Equation  14  shows  the  final 
calculation  of  the  vibrational  energy  where  Nmax  is  the  level  at  which  the  discrete  Boltzmann  distribution  is 
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truncated,  R  is  the  universal  gas  constant,  ^  is  the  ratio  of  vibrational  temperature  to  characteristic  vibrational 
temperature  shown  in  Eq.  16,  and  P,;  is  the  probability  of  a  particle  having  the  ith  level  of  vibrational  energy. 
Assuming  vibrational  energy  is  modeled  as  a  simple  harmonic  oscillator,  the  probability  of  a  particle  being  in  a 
vibrational  level  can  be  calculated  using  Eq.  15.  The  maximum  level,  Nmax,  is  chosen  such  that  the  probability 
of  a  particle  having  a  vibrational  energy  greater  than  that  level  is  less  than  1  x  10-8. 


Nma.x 

Evm  =  J2  Pii6vmR 

(14) 

i= 0 

Pi  =  exp  [-i?]  (1  -  exp  [-?]) 

(15) 

Tv  IB 

*  =  a 

UVIB 

(16) 

At  higher  vibrational  temperatures,  such  that  ?  >  0.2,  vibrational  energies  can  be  sampled  from  the  discrete 
Boltzmann  probability  density  function  without  adversely  affecting  the  efficiency  of  the  MPC  method  and  may 
be  necessary  for  physical  processes  that  directly  depend  on  the  vibrational  energy  distribution  function  such  as 
chemistry.  Based  on  the  results  shown  in  Fig.  21,  a  switching  parameter,  <;,  may  be  used  to  change  from  assigning 
average  energies  at  low  temperatures  to  sampling  energies  directly  from  the  discrete  Boltzmann  energy  probability 
density  function  at  higher  temperatures.  The  current  results  shown  in  Fig.  21  suggest  that  a  switching  value  of 
0.5  may  be  sufficient,  such  that  the  average  value  is  assigned  when  c  <  0.5  and  particle  vibrational  energies  are 
sampled  from  the  Boltzmann  probability  density  function  when  >  0.5.  Regardless  of  the  method  used  to  assign 
vibrational  energies,  discrete  vibrational  levels  are  selected  after  vibrationally  inelastic  collisions. 

An  example  of  results  obtained  by  the  MPC  method  for  Mach  20  flow  over  a  sting-mounted  planetary  probe 
with  a  global  Knudsen  number  of  0.01  is  shown  in  Figs.  22-27.  This  flow  condition  corresponds  to  a  low  density 
experiment  performed  at  the  CNRS  facility  in  France  [58,  59,  60].  Due  to  the  very  low  free  stream  densities 
associated  with  this  expansion  tunnel,  the  vibrational  energy  is  assumed  to  be  frozen  at  room  temperature  for 
all  flow  conditions.  The  variation  in  mean  free  path  around  the  probe  and  initial  and  final  rarefied-continuum 
interface  locations  are  shown  in  Fig.  22.  Two  MPC  simulations  were  applied  to  this  flow  condition,  one  with 
rotational  nonequilibrium  only  applied  in  the  DSMC  method  (Rot.  Eq.)  [28]  and  one  simulation  with  rotational 
nonequilibrium  modeled  through  the  entire  flow  region  (Rot.  Neq.)  [29].  Both  MPC  predictions  are  compared  to 
full  CFD  and  full  DSMC  simulation  results.  Figure  23  compares  the  rotational  temperature  contours  predicted 
by  the  MPC  (Rot  Neq),  full  DSMC,  and  full  CFD.  In  general,  the  MPC  method  has  improved  agreement  with 
full  DSMC  results  over  the  initial  CFD  result  across  the  entire  flow  field.  Even  in  regions  where  the  CFD  module 
is  used,  the  MPC  method  has  obtained  improved  boundary  conditions  from  the  DSMC  solver  which  has  allowed 
the  CFD  solver  to  shift  its  result  to  provide  excellent  agreement  with  full  DSMC. 

Figure  24  illustrates  the  prediction  of  flow  field  properties  made  by  full  DSMC,  full  CFD,  and  both  MPC 
simulations  along  the  extraction  lines  shown  in  Fig.  22.  Vertical  lines  denote  the  interface  location  for  the 
corresponding  MPC  result.  Although  both  MPC  results  remain  in  excellent  agreement  with  DSMC  very  near  the 
surface  along  C2,  the  MPC  results  with  a  single  temperature  modeled  within  the  CFD  module  can  not  accurately 
model  the  strong  thermal  nonequilibrium  that  exists  along  nearly  the  entire  extraction  line.  In  contrast,  the 
MPC  prediction  with  the  ability  of  modeling  a  separate  rotational  temperature  within  the  CFD  solver  remains 
in  excellent  agreement  with  full  DSMC  along  the  entire  extraction  line.  In  addition,  the  DSMC  region  simulated 
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using  the  MPC  method  with  the  ability  to  model  rotational  nonequilibrium  within  the  CFD  solver  is  smaller  due 
to  the  less  restrictive  breakdown  parameter.  Even  though  there  is  an  increase  in  computational  cost  associated 
with  solving  an  additional  conserved  variable  in  the  continuum  solver,  the  reduction  in  cost  of  the  DSMC  module 
outweighs  this  increase  resulting  in  a  decrease  in  the  total  computational  cost  of  the  MPC  method  relative  to 
full  DSMC.  This  has  been  seen  for  other  flows  as  well  [38]. 

Figures  25  and  26,  respectively,  compare  velocity  and  rotational  energy  probability  density  functions  predicted 
by  full  DSMC,  full  CFD,  and  the  MPC  method  at  point  A  shown  in  Fig.  23.  Due  to  the  high  degree  of  collisional 
nonequilibrium  within  the  shock,  the  CFD  velocity  distribution  function  which  is  generated  from  gradients  and  the 
first  order  Chapman-Enskog  expansion  does  not  contain  sufficient  information  to  correctly  generate  the  velocity 
distribution  function  predicted  by  full  DSMC.  In  contrast,  the  MPC  method  is  able  to  remain  in  very  good 
agreement  with  DSMC  throughout  the  entire  velocity  space.  Despite  the  macroscopic  rotational  temperature 
predicted  by  CFD  being  within  5%  of  the  full  DSMC  result,  the  rotational  energy  distribution  function  predicted 
by  full  CFD  is  in  poor  agreement  with  the  DSMC  result  throughout  the  entire  rotational  energy  space.  Similar 
to  the  velocity  distribution  function,  the  MPC  method  remains  in  excellent  agreement  with  full  DSMC  results 
for  the  rotational  energy  probability  density  function. 

Figure  27  compares  the  surface  coefficient  of  heat  flux,  defined  in  Eq.  17,  predicted  by  full  DSMC,  full  CFD,  and 
the  MPC  method  with  available  experimental  measurements. 


Ch 


(17) 


Along  the  fore  body  where  the  flow  is  highly  collisional,  all  three  methods  are  in  good  agreement  with  each  other 
and  the  experimental  measurements.  Despite  this  highly  collisional  flow,  CFD  still  slightly  over  predicts  DSMC 
and  MPC  results.  This  is  due  to  the  inability  to  correctly  model  the  Knudsen  layer  within  the  CFD  solver.  As 
the  flow  expands  around  the  corner,  full  CFD  over  predicts  both  DSMC  and  experimental  measurements  by  over 
an  order  of  magnitude.  In  contrast,  the  MPC  method  remains  in  good  agreement  with  both  the  experimental 
measurements  and  the  full  DSMC  results.  Similarly,  the  MPC  method  remains  in  excellent  agreement  with  full 
DSMC  and  experimental  measurements  along  the  sting  mount,  while  full  CFD  over  predicts  full  DSMC  along 
the  entire  sting  mount. 


Abbate  et  al.  [61]  have  proposed  a  hybrid  method  very  similar  to  that  proposed  by  both  Wu  et  al.  [50]  and 
Schwartzentruber  and  Boyd  [54],  but  with  the  capability  of  simulating  unsteady  flows.  It  has  been  found  that 
the  unsteady  hybrid  method  results  are  in  very  good  agreement  with  full  DSMC  simulation  results  for  unsteady 
shock  tube  problems,  while  requiring  less  than  20%  of  the  full  DSMC  computational  time.  This  hybrid  method 
has  also  been  applied  to  unsteady  expansions  into  a  very  low  density  tank. 

Recently,  Kessler  et  al.  [62]  have  proposed  a  Coupled  Multiscale  Multiphysics  Method  (CM3)  for  simulation 
of  internal  micro  scale  flows.  Here,  the  DSMC  method  is  used  as  a  corrector  to  the  momentum  and  energy 
transfer  terms  in  rarefied  regions.  The  continuum  Navier-Stokes  equations  are  solved  throughout  the  flow  field, 
then  in  regions  identified  as  being  rarefied,  DSMC  particles  are  created  and  assigned  velocities  consistent  with 
the  Chapman-Enskog  velocity  distribution  function.  Next,  the  DSMC  particles  are  allowed  to  move  and  collide 
using  standard  techniques  until  the  initial  distribution  function  has  evolved.  Finally,  the  shear  stress  tensor 
and  heat  flux  vectors  are  updated  in  regions  where  DSMC  data  is  available  and  used  in  the  continuum  module. 
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Figure  22:  Interface  location  and  variation  of  mean  free  path  around  a  hypersonic  probe  [29] 


This  process  is  repeated  until  the  full  solution  has  reached  steady  state.  The  CM3  method  has  been  applied  to 
Couette  flow  over  a  range  of  Knudsen  numbers.  For  example,  Fig.  28  shows  a  schematic  of  the  geometry  used 
to  test  the  CM3  method,  while  Fig.  29  shows  an  example  of  the  simulation  results  for  a  Rayleigh  flow  velocity 
profile  with  a  global  Knudsen  number  of  0.2  based  on  the  channel  height.  This  flow  condition  is  near  the  onset 
of  the  transition  regime  where  the  Navier-Stokes  equations  with  no-slip  boundary  conditions  are  not  adequate 
for  modeling  the  flow,  but  solutions  obtained  using  slip  boundary  conditions  do  produce  reasonably  accurate 
solutions.  With  the  exception  of  the  very  near  wall  region,  the  CM3  method  has  been  successful  at  matching 
DSMC  and  Navier-Stokes  with  slip  predictions.  Improvement  of  the  physical  accuracy  of  the  CM3  near  walls 
and  increased  efficiency  are  ongoing  research  topics. 


5  Remaining  Challenges  and  Summary 

Despite  the  sophisticated  capabilities  of  hybrid  particle-continuum  simulation  techniques,  additional  advance¬ 
ments  of  the  technology  can  still  be  achieved.  One  area  is  further  study  of  a  predictor  of  continuum  breakdown. 
Many  of  the  switching  parameter  values  outlined  in  Sec.  3.1  have  been  found  by  comparison  of  predictions  of 
fully  continuum  and  kinetic  simulation  techniques.  However,  in  all  flows  of  interest  for  application  of  hybrid 
particle-continuum  methods,  the  continuum  regions  are  dependent  on  accurate  prediction  of  nearby  rarefied  re¬ 
gions.  Therefore,  full  simulations  of  these  flows  with  a  continuum  method  can  shift  the  results  in  regions  that  are 
truly  continuum.  Additional  study  of  the  rarefied-continuum  breakdown  and  switching  parameter  using  a  hybrid 
particle-continuum  simulation  technique  may  result  in  a  more  appropriate  switching  parameter  that  reduces  the 
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Figure  23:  Comparison  of  translational  temperature  by  DSMC,  CFD,  and  the  MPC  method  [29] 


computational  cost  of  hybrid  techniques  while  maintaining  the  same  level  of  physical  accuracy. 

Further  development  of  the  hybrid  methods  to  simulate  full  three-dimensional  flows  may  have  the  largest  compu¬ 
tational  savings  over  fully  kinetic  methods.  However,  very  few  hybrid  particle-continuum  methods  have  demon¬ 
strated  the  capability  of  simulating  these  flows  over  complicated  three-dimensional  geometries.  Further  work  in 
this  area  to  mature  the  technology  is  still  required. 

As  seen  with  the  inclusion  of  rotational  and  vibrational  energy  nonequilibrium  within  the  MPC  method,  the 
use  of  additional  physical  models  within  a  hybrid  particle-continuum  code  may  require  extra  considerations  to 
maintain  a  high  level  of  physical  accuracy.  Very  few  studies  of  the  effect  of  dissociation  [63,  64]  and  ionization  on 
continuum  breakdown  have  been  performed.  Among  other  unknown  difficulties,  the  mitigation  of  the  statistical 
scatter  associated  with  trace  species  may  be  necessary. 

Hybrid  particle-continuum  simulation  techniques  have  evolved  over  the  past  two  decades  as  a  powerful  analysis 
tool  for  computation  of  flows  that  can  not  be  fully  simulated  with  either  continuum  or  kinetic  methods  and  still 
maintain  both  physical  accuracy  and  numerical  efficiency.  Typical  multi-scale  flows  that  may  require  application 
of  a  hybrid  particle-continuum  method  have  been  described.  A  review  of  demarcation  and  coupling  procedures 
used  to  construct  hybrid  methods  has  been  provided.  A  review  of  proposed  hybrid  particle-continuum  methods 
along  with  examples  of  a  selection  of  hybrid  simulation  results  has  been  given.  Finally,  a  survey  of  some  of  the 
remaining  challenges  have  been  described. 
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Surface 


Figure  24:  Temperature  and  density  predicted  by  DSMC,  CFD,  MPC  (Rot.  Neq.),  and  the  MPC  method  (Rot. 
Eq.)  along  C2  [29] 
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flat  plate  at  a  20°  angle  of  attack  [62] 
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